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A new numerical scheme to solve the Einstein field equations based upon the generalized har- 
monic decomposition of the Ricci tensor is introduced. The source functions driving the wave 
equations that define generalized harmonic coordinates are treated as independent functions, and 
encode the coordinate freedom of solutions. Techniques are discussed to impose particular gauge 
conditions through a specification of the source functions. A 3D, free evolution, finite difference 
code implementing this system of equations with a scalar field matter source is described. The 
second-order-in-space-and-time partial differential equations are discretized directly without the use 
first order auxiliary terms, limiting the number of independent functions to fifteen — ten metric 
quantities, four source functions and the scalar field. This also limits the number of constraint 
equations, which can only be enforced to within truncation error in a numerical free evolution, 
to four. The coordinate system is compactified to spatial infinity in order to impose physically 
motivated, constraint-preserving outer boundary conditions. A variant of the Cartoon method for 
efficiently simulating axisymmetric spacetimes with a Cartesian code is described that does not use 
interpolation, and is easier to incorporate into existing adaptive mesh refinement packages. Pre- 
liminary test simulations of vacuum black hole evolution and black hole formation via scalar field 
collapse are described, suggesting that this method may be useful for studying many spacetimes of 
interest. 



I. INTRODUCTION 

One of the primary goals of numerical relativity today 
is to solve for astrophysical spacetimes that are expected 
to be strong sources of gravitational wave emission in the 
frequency bands relevant to current and planned gravi- 
tational wave detectors. Expected sources include the 
inspiral and merger of compact objects, supernovae, pul- 
sars, and the big bang. An important tool for extracting 
physics from detector signals is the technique of matched 
filtering, which requires an accurate wave form of a model 
of the expected source. For binary black hole mergers (in 
particular) it is thought that numerical relativity is the 
only method that will be able to provide such waveforms 
close to and during the plunge phase of the merger. De- 
spite significant progress made over the past decade, a 
full solution to this problem still eludes researches. One 
reason for the difficulty is the complexity of the field 
equations. This translates into significant computer re- 
sources being needed to solve the equations, which limits 
the turn-around time for testing new ideas. However, 
perhaps the largest obstacle so far has been finding a for- 
malism to write the field equations in that is amenable 
to long-term, stable numerical evolution. Some of the 
promising techniques used today include symmetric hy- 
perbolic formalisms Q, 0, the BSSN formalism (some- 
times referred to as the NOK formalism) 0, 0, IE El an d 
characteristic evolution (for black hole/neutron star sys- 
tems) 0- Several groups are also beginning to examine 
the possibility of constrained evolution for the 3D bi- 
nary black hole problem 0, 0, 0] , and other promising 
directions make use of tetrad formulations of the field 
equations E^EJEj, and solution of the conformal field 
equations [P^llrllrA ITsj. 

A method of writing the field equations that has proven 
very useful in analytic studies is arrived at by impos- 



ing the harmonic coordinate condition, where the four 
spacetime coordinates x^ are chosen to individually sat- 
isfy wave equations: \3x^ — 0. The Einstein equations, 
when written with this condition imposed, take on a 
mathematically appealing form where the principal part 
of each partial differential equation satisfied by a metric 
component g a p becomes the scalar wave operator Ug a p. 
This allowed for (among other things) the first existence 
and uniqueness proof of solutions to the field equations 
[T^|. In numerical relativity, a solution scheme based 
directly upon this formulation of the field equations has 
recently been suggested by Garfinkle [2(j (see also related 
work b y S zilagyi and Winicour [2l| , and the so-called Z4 
svstem[22j. which seems to be quite similar to generalized 
harmonic evolution in many respects). Garfinkle consid- 
ered a generalization of the harmonic coordinate condi- 
tion of the form Ox* 1 = H^, where are now arbitrary 
source functions, and found that the technique was suc- 
cessful in simulations of the approach to the singularity 
in certain cosmological spacetimes. 

One purpose of this paper is to begin to investigate the 
use of the generalized harmonic decomposition in asymp- 
totically flat spacetimes. The formalism is described in 
Sec. |n] If this method is to be useful for a large class of 
spacetimes, one issue that needs to be addressed is how 
to choose gauge conditions via specification of the source 
functions H^; this topic is discussed in Sec. IIIII A second 
goal of this paper is to investigate direct discretization 
of the second-order-in-space-and-time partial differential 
equations 1 (in other words, the system is not converted 



Recent analytic investigations by Calabrese l23l have suggested 
that such a scheme may suffer from high-frequency instabilities 
in situations where the coefficients in front of mixed time-space 
derivatives are greater than the local characteristic speed. We 
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to a system of first order equations before discretization) . 
One reason for doing so is to have a free evolution scheme 
where the only constraints amongst the variables are the 
four constraint equations imposed by the Einstein equa- 
tions (see also 0, |24|)- The hope then is that even 
if this system suffers from "constraint violating modes" 
2 , it may be easier to analyze and cure them using (for 
instance) ideas suggested in research of symmetric hyper- 
bolic versions of the field equations pi Hi I27L I2I l^ 3 . 
The numerical code is described in Sec. IIVI along with 
related topics such as apparent horizon finding, excision, 
boundary conditions, initial conditions, and the current 
scalar field matter source. Also described in Sec. IIVI is 
a variant of the Cartoon method §£2| to efficiently simu- 
late axisymmctric spacetimes with a Cartesian code. The 
advantages of the method presented here are that no in- 
terpolation is used, and the axisymmctric simulation is 
performed on a two-dimensional slice of the Cartesian 
grid. In Sec. EJtest simulations of black hole evolution 
and gravitational collapse are shown, suggesting that this 
solution method holds promise for simulating asymptot- 
ically flat spacetimes. Concluding remarks are given in 
Sec. IVII in particular a discussion of some of the work 
that still needs to be done before the code could provide 
new physical results in situations of interest. 



II. THE EINSTEIN FIELD EQUATIONS IN THE 
GENERALIZED HARMONIC DECOMPOSITION 



Consider the Einstein field equations in the form 

R aB = 4tt (2T a0 - g a0 T) , (1) 

where R a0 is the Ricci tensor, g a0 is the metric tensor, 
T a0 is the stress energy tensor with trace T, and units 
have been chosen so that Newton's constant G and the 
speed of light c are equal to 1 . The Ricci tensor is defined 
in terms of the Christoffel symbols T^g by 
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where T 1 a is 
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(3) 



have not yet noticed such an instability, probably because of the 
numerical dissipation we us e, which was one of the suggested 
cures for the problem in l23l . 

2 By constraint violating mode we mean a solution to the contin- 
uum evolution equations that is not a solution of the full Einstein 
equations, and furthermore exhibits exponential growth from ini- 
tial data with arbitrarily small deviations from putative initial 
data that does satisfy the constraints. 

3 For it appears that it may not be possible to construct a 
constrained-transport type numerical evolution scheme that sat- 
isfies all of the Einstein equations to machine precision l.'jil l.'ill . 



The notation f iCt and d a f is used interchangeably to de- 
note ordinary differentiation of some quantity / with re- 
spect to the coordinate x a . 

Introduce a set of four source functions via 



H' 
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(4) 
(5) 

(6) 



or equivalently, defining — g fJlU H 1 ' ', we have 

Hp = (hi V^^-gTguw (?) 
The symmetrized gradient of is thus 

H( P ,„) = (In V^g),^ v ~ 9 al3 (,u9^p, a ~ 9 a0 9i3(n,v) a (8) 

The generalized harmonic decomposition involves replac- 
ing particular combinations of first and second deriva- 
tives of the metric in the Ricci tensor (|2J) by the equiva- 
lent quantities in (|7I8|) . and then promoting the source 
functions to the status of independent quantities. 
Specifically, one can rewrite the field equations as 

5* 7 5a/3,7<5 + g l5 ,p9aS,j + 9 1& \a9/3S..j + 2 #(a,/3) 

-2H s T s ad + 2Tj T s ia = -Stt (2T a0 - g a0 T) (9) 

As are now four independent functions, one needs 
to provide four additional, independent differential equa- 
tions to solve for them, which we write schematically as 



^pHp = (no summation). 



(10) 



is a differential operator that in general can depen- 
dent upon the spacetime coordinates, the metric and its 
derivatives, and the source functions and their deriva- 
tives. Note however that the principal part of (JHJ) is now 
the simple wave operator g 5l d 1 ds acting upon each met- 
ric component g a0 \ this subsystem of equations is man- 
ifestly hyperbolic given certain reasonable conditions on 
the metric 4 and as long as the coupling between JJjJ , i|l(J|) 
and any matter evolution equations that may be needed 
do not the affect the characteristic structure of ^fy. We 
will not discuss the well-posedness of this system of equa- 
tions here, though this is certainly a topic worth pursu- 
ing. 

The Einstein field equations are thus equivalent to the 
system of equations © and (fTU|) , provided that the har- 
monic "constraints" (JJJ) are satisfied for all time t = x°. 
The claim then is, at the analytical level, if © is used to 
evolve g a g, and i|10[) is used to evolve H tll then @ will 



4 For example one would need a single coordinate to be timelike 
and the rest to be spacelike throughout the integration volume. 
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be satisfied for all time provided that initial conditions 
are specified so that (Q and © are satisfied then. For 
the special case where are given as a-priori functions 
of the coordinates x^, the preceding statement has been 
proven before [33j (the case = was first shown in 
jlflj)- The idea behind the proof is as follows. Define the 
harmonic constraint function C M as 

C = - Ox". (11) 

For any solution to the Einstein equations , C M is iden- 
tically zero. Using the contracted Bianchi identity and 
conservation of stress energy, one can show that sat- 
isfies the following homogeneous wave equation 

□C = -R^ V C V . (12) 

Therefore, given any that satisfies (0 for all time 
together with some if M that satisfies both = and 
(9 t C M = at t = 0, l|12|) guarantees that g M „ will also 
solve the Einstein equations (JTJ for all time. We can- 
not prove such a result for a general evolution system 
where £f M is specified via some arbitrary set of differen- 
tial equations. Rather, we will take the more pragmatic 
approach in the numerical code of demonstrating con- 
vergence to a solution of the Einstein equations for any 
particular evolution system we use. In fact, such a con- 
vergence test is the only measure of the validity of the 
numerical solution, regardless of any analytic properties 
of the underlying continuum problem. 

Equivalent to enforcing © at t = is to make sure 
that the usual constraint equations, namely 

(3) i? + K 2 - K ab K ab = 167rp, (13) 
K a " \b-K\a = ^J a (14) 

are satisfied then, which from a practical standpoint is 
easier to solve than (JHJ) using existing, well-established 
techniques 0, ^2- In the above, K a b is the extrinsic cur- 
vature tensor of the t = const, hypersurface with induced 
metric h a b, K is the trace of K a t, ( 3 >R is the Ricci scalar 
of h a b, | denotes the covariant derivative operator com- 
patible with h a b, and p and J a are the projected matter 
energy and momentum densities respectively. Note that 
we use notation where Greek indices denote four dimen- 
sional quantities and run from to 3, and Latin indices 
denote three dimensional spatial quantities and run from 
1 to 3. 

At this stage the system (|9I10|) is completely general in 
that we have not yet specified any time slicing or spatial 
coordinates. Choosing a gauge amounts to specifying a 
set of source functions through i|10|) , and thus the source 
functions play a role analogous to the lapse function and 
shift vector in the traditional ADM decomposition. One 
disadvantage of the harmonic decomposition is that (to 
my knowledge) there is no simple geometric description of 
the relationship between and the resulting spacetime 
coordinates. In same cases one can appeal to the ADM 
lapse and shift view of coordinate freedom to motivate a 



particular choice of H^. We will discuss these and several 
other classes of gauge conditions that may be useful for 
numerical evolution in the following section. 



III. SPECIFYING A GAUGE 

Within the generalized harmonic decomposition one 
can think of the source functions as representing the 
four coordinate degrees of freedom available in general 
relativity. There are many conceivable ways of choosing 
H^; in this section we will give a few suggestions, several 
of which are used in the evolutions presented in Sec. W\ 
However, the discussion here is rather heuristic in that 
we do not consider how any of these gauge choices may 
affect the character of the coupled Einstein-gauge evolu- 
tion system. Note that gauge source functions were dis- 
cussed by Friedrich^3| in some detail, though not specifi- 
cally within the context of supplying additional evolution 
equations for them. 

The simplest gauge choice in this formalism is to set 
the source functions equal to some arbitrary functions of 
the spacetime coordinates: 

H M = U(x a ). (15) 

The case / M = is standard harmonic coordinates. The 
next condition we consider is a coordinate system that 
evolves toward harmonic coordinates: 

dtH^ = —K^{t)H^ (no summation), (16) 

where k m are a set of 4 arbitrary though positive functions 
of time, which if non-zero will cause to evolve to zero. 

A useful method to derive coordinate conditions for the 
harmonic decomposition is to appeal to the manner in 
which the coordinate system is specified in the ADM de- 
composition. This makes available a tremendous amount 
of research that has gone into gauge related issues for 
ADM-based evolution @. In the ADM formalism the 
metric element is written as 

ds 2 = -a 2 + hij (dx i + I3 l dt) (dx 3 + ftdt) , (17) 

where the lapse function a measures the rate of change 
of proper time with respect to coordinate time t of hy- 
persurface normal observers, is the intrinsic metric of 
t = const, slices, and the shift vector (3 J describes how 
the spatial coordinates change for normal observers from 
one time slice to the next. The normal component and 
spatial projection of the source function are 5 

H-n = 



5 Note however that does not transform like a one-form under 
coordinate transformations, and hence the projections are not 
covariant objects. 
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-K-d u Qna)n v 



(18) 



= -T*j k h? k + djQxia)^ + -d^ffn 1 , (19) 
J a 

where T l - k is the connection associated with hij, and n v 
is the hypersurface normal vector given by 



= —ad v t 



(20) 



Notice that in l|18f) and l|19f) the time derivative of a only 
appears in H ■ n, while the time derivative of /3 Z only ap- 
pears in the corresponding component of _L H l . In other 
words, the choice of the normal component H-n in an evo- 
lution directly affects the rate of change of a with respect 
to time, and therefore H-n controls the time-slicing of the 
spacetime; similarly, _L H l controls the manner in which 
the spatial coordinates evolve with time (another way of 
stating this is that (|18I19|1 are generalizations of the hy- 
perbolic equations governing the lapse and shift within 
harmonic coordinates ^ij). One way in which an ADM 
style gauge condition can be used within the harmonic 
decomposition is to substitute the corresponding choices 
of a and (3 % into (|18I19|I , and use the result as the source 
functions for the harmonic evolution. The simplest class 
of gauge conditions that can be implement in this fashion 
are the so called "driver" conditions [iBlliR I4 H II I |49 | . 
where one directly specifies the time derivatives of a and 
[3 l to achieve, for example, approximate maximal slicing 
and minimal distortion gauges respectively. 

The manner in which the ADM driver conditions are 
implemented suggests a similar way in which such gauge 
conditions can be used in a harmonic evolution: instead 
of substituting in the forms for a and [3 l in i|18ll9fl to try 
to satisfy the conditions exactly, choose source functions 
to drive the gauge toward the desired one. To see how this 
can be done, first rewrite H18I19|I as evolution equations 
for the the lapse and shift: 



dtct — —a 2 H ■ n 
dt? = a 2 _L HU 



(21) 
(22) 



where the ellipses denote the rest of the terms that do not 
contain dta, dtff or H^. Now suppose at some instant of 
time the desired value of the lapse and shift are calculated 
(by whatever means) to be ao and (3q respectively. Then 
from (|21I22() one possible set of choices for the source 
functions that will cause the lapse and shift to evolve 
toward the desired values are 



H-n 
_L H l 



K n (t) 



a — ao 



-Ki(t) 



a 



(23) 

(no summation), (24) 



where n n and are positive functions of time that can 
be used to control the rate of evolution. 

It many circumstances it may make more sense to im- 
plement the above style driver conditions as evolution 



equations, rather than algebraic conditions. This could 
be, for example, if the initial conditions for are not 
compatible with the desired gauge choice, and so imple- 
menting (|23I24|) will result in discontinuous source func- 
tions at the initial time. A couple of alternative possibil- 
ities include 



dH ■ n a~ a$ 

at a 1 



(25) 



and 



U{H ■ n) = - Kn (t)^^ + f n (t)(H ■ n) l|t n", (26) 



with similar expressions for the spatial parts of H^ . (t) 
is a positive function that can be used to add a dissipative 
term to (12611 . One advantage of using a wave operator 
(|26|l to evolve the source functions is then the principal 
parts of all equations in the system I|9I1U|I have the same 
characteristic structure (as long as ao and /3q depend at 
most on first derivatives of the fundamental variables). 
This may be important to establish well-posedness of the 
coupled system of equations [50| . 

It is beyond the scope of this paper to investigate how 
well any of these suggested gauge conditions perform in 
situations of interest, however in Sec. [V] we will show 
some preliminary results indicating that these ideas can 
be implemented in a stable fashion. 



IV. NUMERICAL CODE 

In this section we describe a 3D numerical code based 
upon the generalized harmonic decomposition. This code 
has several features of note 

• Direct discretization of f$ llU\) : In other words, we 
do not convert the system of equations to first order 
form — the only variables used are the 10 unique 
metric components g a /3, 4 source functions H^, and 
matter variables. 

• Spatially compactified Cartesian coordinates: The 
spatial coordinates are compactified to include i°, 
to simplify the imposition of physically realistic 
boundary conditions for asymptotically flat space 
times. 

• Black hole excision: Black hole excision is used to 
evolve spacetimes containing black holes, whereby 
portions of the computational domain inside of ap- 
parent horizons are "excised" to remove the singu- 
larities. 

• Built within a parallel adaptive mesh refinement 
(AMR) framework: The code utilizes a new set of 
parallel AMR libraries which we will describe else- 
where, though the Berger and Oliger style AMR 
algorithm used is very similar to the one presented 
in [Elil. 



• Efficient simulation of axisymmetric spacetimes us- 
ing a variant of the Cartoon method ]33j : The algo- 
rithm presented here does not require interpolation, 
and only utilizes a single 2-dimensional slice of the 
Cartesian grid, simplifying incorporation into ex- 
isting AMR packages. 

In the remainder of this section we will describe certain 
aspects of the code in more detail. 
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A. Discretization Scheme 

From here onward we will use the coordinate names 
t = x a and (x,y,z) = (x 1 , x 2 , x 3 ). Also, as discussed in 
the next section, we use a compactified coordinate sys- 
tem in the code. This necessitates the use of regularized 
variables for some of the metric and source function com- 
ponents, however to keep the discussion in this section 
simpler we ignore that aspect of the code here. In Ap- 
pendix[B]we present a stability analysis of this discretiza- 
tion method applied to a one-dimensional wave equation 
in flat space. The purpose of the analysis is to give a sim- 
ple, concrete example of the numerical method, and to 
show that there are no fundamental instabilities in it. Of 
course, this cannot prove that the full, non-linear prob- 
lem in compactified coordinates will be stable — doing so 
is beyond the scope of this paper. 

We use second order accurate finite difference tech- 
niques to discretized (|9I1U|I and the scalar field evolution 
equation 1|45|) presented in Sec. IIVDI This is a set of 15 
equations for 15 unknown functions — the 10 non-trivial 
metric components g a /3, the 4 source functions and 
the scalar field <&. In the discretized version of @ all 
Christoffel symbols, contravariant metric elements and 
their gradients are replaced with the appropriate sum of 
covariant metric elements and their gradients. As © 
and H45fl are second order partial differential equations 
in time, second order accurate discretization requires a 
3 time level scheme (at a minimum) . Fig. ^ shows a 
schematic representation (with two spatial dimensions 
suppressed) of the discretization of a variable f(t, x, y, z). 
f(t,x,y,z) evaluated at a grid location (t n ,Xi,yj, Zk) = 
(nAt, iAx,jAy, kAz) is denoted by / 8 " fe , where n, i, j and 
k are integers, and At, Ax, Ay and Az are the temporal 
and spatial discretization scales respectively. Table^con- 
tains a representative sample of the finite difference oper- 
ators used to evaluate derivatives on the mesh. Replacing 
the continuum variables with discrete variables, and the 
derivative operators with difference operators will result 
in a difference equation 



C 



flijk 







(27) 



for each variable / at each grid point 
(t, x, y, z) = (t n , Xi, yj, Zk) in the computational domain. 

We solve the system of equations (|27H using a Newton- 
Gauss-Seidel relaxation scheme, as follows. Initial data 
for a single time step at t = t n consists of all the variables 



-H- 



-B — 



FIG. 1: The discretization of a variable f(t,x,y,z) in the x — t 
plane. 



at time levels t n and t n 
ables at time level t n+1 . 



L . The unknowns are the vari- 
Denote an approximate value 



of the unknown f™^ 1 by f^A ■ The iteration is set up 
using function values at time level n as an initial guess to 
the solution at time level n + 1. One step of the iteration 
then proceeds by updating each unknown, in turn, via 



fn+i 
Jijk 



in+l 
Jijk 



K 



f\ijk 



Jf 



(28) 



tjk 



where TZf is the residual of the difference equation 
(the left hand side of l|27|) evaluated using the approxi- 
mate solution) and Jf\™j k is the diagonal element of its 
Jacobian 



3f% 



dC 



ijk 



f\ijk 



df^k 1 



(29) 



again evaluated with the approximate solution. In other 
words, (|28|) is simply solving a linearized version of l|27|) 
for Z™-^ 1 assuming all other unknowns are fixed. The 
iteration is repeated until the residual for all variables is 
below some specified tolerance. 



1. Numerical Dissipation 

Some form of numerical dissipation is necessary to sta- 
bly evolve certain spacetimes, in particular those that 
contain black holes. We use Kreiss-Oliger style dissipa- 
tion |53j| . however, rather than modify the discrete evolu- 
tion equations as is typically done (and note also that |53| 
considered first order in time systems) , we apply the dis- 
sipation as a filter to the discrete variables, at both past 
time levels t n and i n_1 , prior to updating t n+1 . Specifi- 
cally, at a given time level we define the high-frequency 
component of grid function fij k , in the x direction 
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f,t 

f,XX 
f.tt 



(fill ~ /T-i)/(2Ax) 



n+l 



/r~ 1 )/(2At) 



(/ „+l_ 2/ n + / n-l )/(At) 2 



/• ( f n+1 
J ,tx VJj+l 



f, 



n+l 
1 



TABLE I: A sample of the finite difference stencils used to 
convert the differential equations to difference equations. The 
column on the right shows the second order accurate repre- 
sentation (with y and z indices suppressed for clarity) of the 
corresponding derivative operator to the left, evaluated at the 
point (t n ,Xi,yj, Zfc). Similar stencils are used for terms con- 
taining y and z derivatives. 



as 



Vijk 



Y^(/i-2jfe — ^fi-ljk + Sfijk 



~^fi+ljk + fi+2jk), 

elsewhere, 



2 < i < N T - 2 



(30) 



where the local size of the mesh is N x points in the x 
direction. After r]fj k has been calculated over the entire 
local grid, it is subtracted from / as follows 



fijk — fijk e Vijki 



(31) 



where e is a constant, required to be in the range 0..1 
for stability. In practice we use values of e in the range 
0.2 to 0.5. Once the high frequency components in the x 
direction have been subtracted, the procedure is repeated 
for the high frequency components rjf^ and rj^ k of fijk 
in the y and z directions respectively, which are given by 
expressions similar to l|30[l. 

We did experiment with extending the dissipation fil- 
ter to the grid boundaries as outlined in |54j . however 
this did not seem to have a significant effect on the solu- 
tion in most circumstances, and seemed to produce more 
error (as measured by residuals of the field equations) 
next to excision boundaries without offering improved 
stability. However, the excision method proposed in [54| 
was for cubical excision boundaries, and for schemes sat- 
isfying summation by parts, so it is questionable how 
appropriate it is to apply that method here. Also note 
that applying the above filters to both past time levels 
at each evolution step is essential for long term stability. 
We do not know why this is so important; naively one 
would think that only applying the filter to time level t n 
would be sufficient, as the update step does not alter the 
variables at t n , and since t n is copied to t n ~ 1 after each 
update step, both t n and t n_1 are effectively smoothed. 
Also, a simple extension of the analysis in Appendix [B] 
to account for different amounts of dissipation applied to 



each of the past time levels shows that the one dimen- 
sional, flat space wave equation remains stable; hence the 
need to dissipate both time levels is particular to black 
hole spacetimes as far as we can tell. 



B. Coordinate System and Boundary Conditions 

To simplify the imposition of asymptotically flat 
boundary conditions we use the following spatially com- 
pactified coordinate system. First, consider an uncom- 
pactified Cartesian coordinate system of the form 



ds 2 



gudt 2 + 2g ti dx l dt + gijdx l dx^ 



(32) 



Here (x 1 ,^ 2 ,^: 3 ) = (x 7 y,z) runs from — oo to +oo, and 
in the limit where x l ±oo the metric becomes the 
Minkowski metric 



ds 2 = -dt 2 + dx 2 + dy 2 + dz 2 . 
The following coordinate transformation 



(33) 



(34) 



x l = tan(7nz;'72) 
(with i = t) will bring (|32|l into the form 

ds 2 = gudt 2 + 2g t idx l dt + g^dx 1 dx\ (35) 



where 



gu = -sec 2 (7rx72)3 ti , 
■ sec 2 (7ra72) sec 2 (7rx J 72)5y. 



(36) 



Now x l runs from —1 to 1, and spacelike infinity i° cor- 
responds to the surfaces x l — ±1. Note that in this limit 
the compactified (unbarred) metric elements are singular, 
however the uncompactified parts are still well behaved 
and asymptote to their Minkowski values. In the code we 
thus evolve the uncompactified components gtt,gu and 
cjij , analytically substituting the values (|36|l into © prior 
to discretization. Furthermore, in the compactified coor- 
dinate system (|35|l we define the spatial source functions 
Hi to take the form 



Hi^ Hi- ntaxi(irx i /2). 



(37) 



and evolve only the regularized components Hi (for 
note that in compactified Minkowski coordinates Hi = 
7rtan(7ra;72) from Q)). Therefore, the outer bound- 
ary conditions we impose on the regularized metric and 
source functions are 



9tt(t,i°) 
9u(t,i°) 
gu(t,i°) 
9ij(t,i°) 



-1 



1 

o, i + j 
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H t (t,i°) 
Hi{t,i°) 



by 



(38) 



where the notation (t, i ) refers to any one of the six 
boundaries x = ±1, y — ±1 and z = ±1. 

We conclude this section by discussing several concerns 
about evolving the field equations in a coordinate system 
compactified to spatial infinity. It is beyond the scope of 
this paper to analyze these issues in more detail, how- 
ever, we are currently investigating them. However, note 
that a similar compactification scheme was used to model 
black strings in 5 dimensions [3^, with no notable ad- 
verse effects. 

First, the metric, hence equations, are formally sin- 



gular at 



The singular behavior is dealt with us- 



ing regularized variables and enforcing Minkowski space 
boundary conditions at i°, as described above. Never- 
theless, there are terms in the equations that grow like 
l//i 4 at grid locations near the outer boundary, where 
h is the mesh spacing there. Therefore, for the equa- 
tions to remain regular near i° during evolution requires 
that the leading order behavior of the metric and scalar 
field variables always approach there asymptotic values 
sufficiently fast to cancel this divergent behavior (this 
is essentially the same problem one must deal with in 
an axisymmetric code near the axis singularity) . For the 
simple test results presented in Sec. 0the evolution near 
i° is well behaved, however we cannot guarantee that this 
will be the case for all classes of asymptotically flat initial 
data. 

A second issue is the propagation of outgoing waves 
toward i°. The compactification causes the wavelengths 
and speeds to decrease. Thus, for any fixed resolution 
near i°, such waves will eventually be poorly resolved 
on the grid 6 . This could lead to a couple of undesir- 
able effects. First, numerical dissipation will significantly 
decrease the amplitude of the waves, making waveform 
extraction in the outer regions of the domain impracti- 
cal. Second, some portion of the wave will get "reflected" 
back to the interior of the domain, which is not physi- 
cal and may adversely affect the accuracy of the interior 
solution. 



C. Apparent Horizon Finder and Excision 

We use the following flow method to search for sin- 
gle, simply-connected apparent horizons in the spacetime 
(this is the same algorithm used in l&l: see [3^ for a re- 
view of most current methods, and |4Ct l4l| for some re- 
cent work on fast, elliptic-solver based apparent horizon 
finders). Consider the level set function F(r, 9, </>) defined 



= r-R(9,, 



(39) 



where the spherical polar coordinates (r, 6, </>) are defined 
in terms of uncompactified coordinates, relative to some 
center (x ,y ,z ), via 



x = xq + r cos 4> sin 9 
y = ya + T sin <f> sin 9 
z = zq + r cos 9 



(40) 



We want to find the function R(9, </>) such that the hy- 
persurface F = has zero outward null expansion 



9 



Ka;/3< 



a/3 



(41) 



where h a @ is the spatial metric (|17[) and £ a is the outward 
pointing null vector normal to F = const, surfaces: 



y/h s ^dsFd 7 F' 



(42) 



The flow method involves specifying some initial guess for 
R, then evolving the following equation until the magni- 
tude of the norm of evaluated along F = is as close 
to zero as desired: 



dR(0,< 
d\~ 



-e(r,0,0)| r= 



(43) 



where is evaluated along F = 0. This equation is 
parabolic in "time" A, hence dX must be of order (Ax 1 ) 2 
for stability. 

During a typical evolution where a black hole forms via 
scalar field collapse we initialize R(9, (f) = fo, where ro is 
a constant close to though outside 7 of where we expect 
the apparent horizon (AH) to first form, and periodically 
(every tens to hundreds of time steps) search for an AH 
using this initial guess until one is found. For subsequent 
AH searches we use the previously found surface as an 
initial guess for R(9, 4>). If multiple black holes form we 
search for each AH independently. 

Some form of excision is necessary for long term evo- 
lution of spacetimes containing black holes. Excision 
means that one places interior boundaries inside of all 
black holes such that all physical singularities are re- 
moved from the computational domain. This assumes 
that cosmic censorship holds, which further implies that 
a black hole's event horizon will be outside any appar- 
ent horizon, and hence one can use the apparent horizon 
as a guide where to excise. For each black hole, we ex- 



Keeping the waves well resolved with AMR is not a practical 
solution in general, as the outgoing wavetrain one expects from 
a binary inspiral, for example, is volume filling. 



7 The underlying assumption in 14.'H is that © > implies that 
the surface is outside of the apparent horizon, which is not true 
everywhere at early times during a gravitational collapse simu- 
lation. 
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TABLE II: A sample of the finite difference stencils used to 
convert the differential equations to difference equations ad- 
jacent to an excision surface. The column on the right shows 
the second order accurate representation (with y and z indices 
suppressed for clarity) of the corresponding derivative oper- 
ator to the left, evaluated at the point (t n ,n,yj,Zk)- The 
operators shown above are used when the point Xi-i is inside 
the excision surface and the points Xi,Xi+i, ... are outside of 
it. 



cise along an ellipsoid in compactified coordinate space, 
where the shape of the ellipsoid is chosen to match that 
of the apparent horizon as closely as possible along the 
ellipsoid's principal axis (which currently are required to 
lie along the coordinate axis). The size of the ellipsoid 
is typically a bit smaller than that of the AH, to give 
some buffer zone between the excision surface and the 
AH. Any point on the grid inside the ellipsoid is defined 
to be excised, hence the excised region will necessarily 
be a grid-based approximation to the smooth ellipsoidal 
shape (this is often referred to as "lego excision" in the 
literature) . 

In general, boundary conditions need to be applied 
along the excision surface; however, in a free evolution 
(such as described here) where all the characteristics on 
the excision surface are directed inward, no boundary 
conditions should be placed on the field variables. In 
the current version of the code we assume that this is 
true, though we do not explicitly compute any of the 
characteristics. For a finite difference scheme, such a "no 
boundary" boundary condition means that the evolutions 
equations are applied at the excision surface, with cen- 
tered difference operators replaced, as appropriate, by 
forward or backward difference operators so as not to 
reference grid values inside the excised region. See Table 
[H]for samples of the particular stencils we use. Note that 
we define the excision surface to be constant in time, and 
hence only spatial difference stencils need to be modified. 
During evolution, if the excision surface moves such that 
previously excised points (interior points) become "unex- 
cised" , we initialize them via fourth order extrapolation 
from adjacent exterior points at all time levels in the 
grid hierarchy. We cannot a-priori prove that this exci- 
sion method is stable, rather, as discussed in Sec. [H] wc 
will require convergence to a self-consistent solution of 
the field equations as a proof-by-example that the code 
is stable and correct. 



D. Matter Source 

The present matter source modeled in the code is a 
massless scalar field <f>. The corresponding stress-energy 



tensor T^ u is given by 

T llv = 2$, u p, v -g llv *„&~<, (44) 

and the evolution of $ is governed by the wave equation 

□$ = $./ = 0. (45) 

Note that l|44|) differs by a factor of 2 from the convention 
of Hawking and Ellis [55| , which amounts to rescaling $ 
by a factor of 

E. Scalar field Initial Data 

At this stage for scalar field gravitational collapse we 
only consider time-symmetric initial data with a confor- 
mally flat spatial metric. Specifically, at t — the metric 
and its first time derivatives take the following form: 



9tt(t 
9u(t 
9ij(t 

dtg a p(t 



0,x,y,z 
0,x,y,z 
0,x,y,z 
0,x,y,z 
0,x,y,z 



-1 









! * 7^ 3 
x,y,z) 



(46) 



The scalar field is thus the source of all non-trivial geome- 
try at t = 0, and we initialize it as a sum of Gaussian-like 
functions of the following form 

$(t = 0, x, y, z) = ^2f(x,y,z), 

i 

d t $(t = 0,x,y,z) = 0, (47) 



with 



P{x,y,z) 
P l {x,y,z) 



^expH^y^/A*] 



il2\ 



= {[1 



[1 



V 

A 2 



}[x(x) - xl] 2 + 

[y(y)-yo? + 
Mz)-4]*} 1 ' 2 , 



(48) 



where Aj, A 1 , e*, e 1 , e*, 2q, j/q and Zq are constants, and 
x(x),y(y) and z(z) are given by l|34|) . 

In l|46(l . ^>{x,y,z) is solved for using the Hamiltonian 
constraint (|13[) equation, using an adaptive multigrid 
routine as discussed in [5l], |52( ■ Note however that some 
complications do arise when attempting to solve an ellip- 
tic equation in compactified coordinates using multigrid; 
we will briefly discuss these issues in Sec. IIVE II The mo- 
mentum constraints (|14(l are trivially satisfied with the 
above initial conditions. Once the constraints have been 
solved, we initialize the source functions H a {t = 0, x, y, z) 
using iJSJl and Q. 

With a three time level evolution scheme, the past time 
level at t — —At needs to be initialized as well. To 
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obtain second order accurate convergence of the solution 
at late times, the past time level needs to be consistent 
with the initial data to within At 2 . In the code we have 
implemented a couple of methods to achieve this; the first 
is to use a Taylor expansion along with the equations 
of motion, the second is to evolve backward in time to 
t = —At with a smaller time step. The first method 
works as follows 56] . For any one of the evolved grid 
functions /™- fe the past time level n = — 1 is initialized 
to second order accuracy using a Taylor expansion about 
t = 

„ A/ 2 

fvl=fi jk -fMt + fijl^, (49) 

where f^ k is the first time derivative of f(t,x,y,z) at 
t = (from the initial data), and is the second time 
derivative of f(t, x, y, z) at t = 0, evaluated by substitut- 
ing the initial data into the relevant equations of motion 
(|9I45II and solving for d t d t f- For the second method, the 
past time level is only initialized to first order using 

and ft®}., however with a smaller time step At s w At 2 . 
This initial data is then evolved backward in time un- 
til t = —At, and the solution obtained there is used to 
initialize the past time level for the actual evolution. 



1. Multigrid in a Compactified Coordinate System 

Standard geometric multigrid (MG) methods that use 
point-wise relaxation as a smoother (as we do) are only 
efficient when the size of the coefficient functions multi- 
plying each of the principal parts of the elliptic operator 
are of comparable size [52] . This is the not the case near 
i° in our compactified coordinate system. To illustrate, 
consider the form of the spatial Laplacian V 2 using the 
coordinates of Sec. IIV Bl in flatspace 

V 2 - A ( cos 4 (nx/2) + cos 4 (ny/2) + 

cos 4 (W2)j^)+..., (50) 

where the ... denote lower order terms. Notice that near 
any one of the outer boundaries x % — ±1 the correspond- 
ing coefficient of the second derivative term goes to zero. 
We have not solved the issue of multigrid inefficiency in 
this part of the domain, however if we use scalar field ini- 
tial data of compact support in a sufficiently small region 
about [x, y, z) — (0, 0, 0), then we find that the fine-grid 
relaxation performed within MG is adequate in obtain- 
ing a solution of sufficient accuracy near i (for ^ will 
then go like 1 + 0(1/ r), and the initial guess of \& = 1 is 
close enough to the solution that relatively few relaxation 
sweeps are needed). 

A more serious problem is that relaxation using stan- 
dard centered difference approximations for first deriva- 
tives is unstable near i . A partial solution to this prob- 



lem is to use the following 4-way corner-averaged differ- 
ence operator at grid point (i,j,k) (shown here for the 
x derivative; the other first difference operators are sim- 
ilarly modified) 

f.x = (/i+lj + l,fe+l - /i-l,j + l,fc+l) /(8Ax) + 

(fi+i.j+i,k-i - fi-i,j+i,k-i) /(8Ax) + 
(/i+i,j-i,fc+i - fi-i,j-i,k+i) /(8Ax) + 
(/i+i,j-i,fe-i — /i— i,fe— l) /(8Air) 
+0(Ax 2 ). (51) 

In the limit where the mesh spacing goes to zero in the 
vicinity of a; 1 = ±1, even this modification exhibits relax- 
ation instabilities for initial data that is not sufficiently 
compact about r — 0. However, this is not a problem 
for the kinds of physical systems we plan to use the code 
for, as all the interesting dynamics will be confined to a 
small region abour r = 0, and this will be the part of the 
hierarchy with high resolution. 

F. Exact Schwarzschild Black Hole Initial Data 

For some of the tests describe here we use analytic 
initial data from a Schwarzschild black hole solution 
in Painleve-Gullstrand-like (PG) coordinates. The non- 
compactified components of the metric are 

ds 2 = - M - dt 2 + dx 2 + dy 2 + dz 2 

2\ / 2M 

+ _ 3 ^ 2 [xdx + ydy + zdz] dt, (52) 

where r = ^Jx 2 + y 2 + z 2 and M is the mass of the black 
hole. This (with appropriate spatial compactification) 
gives initial data for the metric at t = 0; and is also 
used to evaluate {7J for the initial values of the source 
functions. 



G. Efficient Simulation of Axisymmetric 
Spacetimes 

In this section a variant of the "symmetry without 
symmetry" , or Cartoon method [32| for efficient evolution 
of an axisymmetric spacetime with a Cartesian-based 3- 
dimensional code is described. The advantages to the 
approach presented here are that no interpolation is ever 
performed, the axisymmetric grid structure is a two di- 
mensional slice of the Cartesian grid, rather than a thin 
three dimensional slab, and the method is not specific to 
finite difference based codes, so can readily be applied to 
a spectral code, for instance. Having the grid structure be 
two dimensional is helpful in that it allows easy integra- 
tion of the code with standard adaptive mesh refinement 
(AMR) packages. The reason is as follows. In the origi- 
nal Cartoon algorithm, the third, thin dimension is one 
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finite-difference stencil- width thick. However, most AMR 
algorithms can only refine a given volume of a grid, which 
would increase the width of the slab-dimension on finer 
levels, and thereby reduce the efficiency of the Cartoon 
method. Of coarse the AMR algorithm could be modi- 
fied to deal with such a situation, however by using a two 
dimensional grid structure one avoids this problem alto- 
gether. Note also that the purpose of the algorithm pre- 
sented here is merely to provide an efficient way to simu- 
late axisymmetric spacetimes with a Cartesian code, and 
not to address any issues of axis stability in axisymmetric 
codes, which was one of the orig inal motivations behind 
Cartoon. Some recent work[34| has suggested that the 
standard Cartoon algorithm may not be stable. Here we 
deal with the axis by applying appropriate regularity con- 
ditions and numerical dissipation, which has proven to be 
an effective method for dealing with the axial singular- 
ity in axisymmetric codes [35ll36| (note also that in some 
cases stability in axisymmetric codes can be obtained by 
constructing schemes with a conserved discrete energy, 
using operators that satisfy summation by parts — see for 
example [33) • 

The idea behind our modified Cartoon algorithm is as 
follows. In a c?-dimensional axisymmetric spacetime we 
have an azimuthal killing vector £ M , hence the metric 
and scalar field matter source $ satisfy 



0. 
= 0. 



(53) 



What these equations imply is that all non-trivial struc- 
ture of the metric and scalar field are encoded within a 
d — 1 dimensional sub-manifold S of the spacetime, as 
long as £ M is nowhere tangent to S. Therefore, one only 
needs to solve the field equations on S, and (|53[1 can be 
used to extend the solution throughout the spacetime. 
In a numerical evolution, it makes most sense to have S 
coincide with a constant coordinate hypersurface, which 
we set to z — for concreteness. We then choose coor- 
dinates such that £ M has the following explicit form (in 
uncompactified coordinates) 



d 



e = y -5= -2 -5= 



dz 



d 



dy 



(54) 



which implies that £ M is orthogonal to z = 0, and the 
axis of symmetry runs along the x direction and is cen- 
tered at y = 8 . To solve the field equations on z — 
requires first and second derivatives of metric variables 
both within the hypersurface 2 = 0, and orthogonal to it 
in the z direction. To calculate z derivatives the original 
Cartoon method effectively extends the solution using 



In other words, 1541 is merely the Cartesian form of (<9/<9</>) M , 
where is a standard azimuthal coordinate with the axis of sym- 
metry coincident with the x axis. 



(|53|l to a sufficient number of grid points above and be- 
low z — 0, so that the usual finite difference stencils can 
be used to calculate z derivatives. The approach taken 
here is to substitute the explicit form of the Killing vector 
(|54|l into the definition l|53l) , and use the resulting expres- 
sion to evaluate the z gradients directly. In other words, 
the same numerical method is used to solve equations 
as outlined in Sec. IIV Al however instead of calculating 
z derivatives using finite different approximation, the z 
derivatives are replace with appropriate combinations of 
x and y gradients. In Appendix lAl we list the results of 
this calculation for all relevant gradients of the metric in 
compactified coordinates; here we illustrate the technique 
for the simpler case of the scalar field <f>. 

Evaluating l|53[) for $ using l|54|) . we obtain 



dz 



z9$ 
V dy 



(55) 



Taking the z derivative of this equation, and replacing 
any z gradients of <E> appearing on the right hand side 
with lJ55)l. gives 



<9 2 $ 

dz 2 



1 <9$ 

y dy 



y 



1 d<S> 

V dy 



<9 2 $ 
dy 2 



(56) 



Evaluating these equations at z = gives 
3$, 



dz 
<9 2 $, 



-o = 0, 

19$ 

dz 2 ~ y dy 
All other mixed second derivatives of <I> involving z are 



(57) 



One thing to note from equation (|57|l is that the axis 
y = is singular. Therefore a regularity condition needs 
to be applied there, which can easily be seen from l|57|) 
to be d$/dy = at y = 0. 



V. PRELIMINARY RESULTS 

In this section we present results from several test sim- 
ulations demonstrating certain aspects of the code. Sig- 
nificantly more work needs to be done before the code 
may be able to produce new physical results, however 
the current simulations suggest that the generalized har- 
monic decomposition could be a viable alternative to the 
ADM decomposition for many problems of interest. 

In Sec. IV Al we show a convergence test of scalar field 
evolution in 3D, Sec. IV Bl evolves a Schwarzschild black 
hole in Painleve-Gullstrand coordinates, and Sec. IV CI 
demonstrates gravitational collapse of scalar field initial 
data to a Schwarzschild black hole. 
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A. Convergence Test 

For a 3D convergence test we used the following ini- 
tial conditions for <!> (|47|) . which describes three prolate 
spheroids slightly offset from one another so that there 
is no spatial symmetry in the problem: 



.4 1 
A 1 

/_2 -2 -2^ 



0.034, A 2 = 0.033, A 3 
0.1 A 2 = 0.1 A 3 = 0.1 
(0.025,0,0), 
(0,-0.025,-0.025), 
(-0.025.025,0.025) 



0.033 



= 0.1, 



< = 0.1, 



= 0.1 



(58) 



and all other initial data parameters for $ are zero. 
With these parameters the ADM mass of the spacetime 
is roughly 0.005, so the initial distribution of energy is 
concentrated in a radius about 10 times larger than its ef- 
fective Schwarzschild radius. For the H t coordinate con- 
dition we used a slightly modified version of H26|) . where 
we eliminated the coupling (through the normal n^) to 
Hi and added an arbitrary power n of a in the denomi- 
nator: 



Off* 



■Kt(t)^nr + 6(^,/. (59) 



where n t (t) = n q{t), £ t (t) = £o?(i), and q(t) is given by 



15- 
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0, elsewhere. 



< t < h 
(60) 



q(t) provides a smooth (twice differentiable) transition 
from at t = to 1 at t = ti, and makes the evolution of 
the source functions consistent with the choice of time- 
symmetric initial data. We evolve Hi to zero using a ver- 
sion of (|16|l with Hi replaced by Hi, and Ki(t) = Koq(t). 
For this particular simulation we had ko = 50, £o = 10, 
n = 3, t\ = 1/10 and ao = 1. 

A convergence test involves running a given simulation 
at several different resolutions, and comparing the results 
to ensure that the solution of the finite difference equa- 
tions is converging to a solution of the partial differential 
equations. We ran three simulations of differing resolu- 
tion. The coarsest resolution run had a base grid size of 
17 3 , and we specified a value for the maximum desired 
truncation error so that up to 5 additional levels of 2 : 1 
refinement were used, giving an effective finest resolution 
of 513 3 — see Fig. below for a depiction of the mesh 
structure at two times during the simulation. We used 
a Courant factor of 0.25 at each level in the hierarchy 
(i.e. At = O^Aa; 1 ). For the medium and finest reso- 
lution simulations we used the same grid hierarchy pro- 
duced by the coarsest resolution simulation (which was 
produced using standard truncation error estimate meth- 
ods), though doubled and quadrupled the resolution of all 





FIG. 2: A depiction of the adaptive mesh structure for the 
convergence test simulation described in Sec. IV Al The image 
to the left corresponds to the mesh structure at t = 0, while 
that to the right at t = 0.5. The largest box in each figure, 
whose faces are at i° , actually represent two levels of (2:1) 
refinement. The increase in size of the finer levels and loss of 
the finest level of refinement by t — 0.5 is due to the outward 
propagation of the initial distributions of energy. 



the grids respectively, keeping the same Courant factor. 
To keep the computational cost of the highest resolution 
run manageable, we only ran the simulation until t = 0.5; 
however this corresponds to roughly five light-crossing 
times of the central region of the grid where the scalar 
field is concentrated, and so a reasonable amount of dy- 
namics does occur. Also, this run time is sufficiently long 
that possible adverse effects from the AMR algorithm, 
such as from regridding or high-frequency "noise" from 
parent-child refinement boundaries, can be captured by 
the convergence test. 

Label some grid function / from the finest resolution 
simulation fh, from the medium one fih and from the 
coarsest one f^. Then the convergence factor Qf we 
calculate is 



Q f 



1 

fn~2 



(^\\f 4h -f 2h \\-ln\\f2h-f h \\), 



(61) 



where before the subtraction we interpolate the grid func- 
tions to a common uniform grid, and then compute the 
£2 norm of the differences. For an n th order accurate 
scheme one would expect Qf to approach a value of n 
in the limit as the mesh spacing goes to zero. See Fig. 
13 for the convergence factors from the above simulations 
for several representative functions. The plot shows that 
we do see convergence close to second order. At early 
times, the convergence factor is slightly worse than sec- 
ond order; we surmise that the reason for this is a small 
amount of unphysical, high-frequency solution compo- 
nents ("noise") present at parent-child mesh refinement 
boundaries at the initial time. This noise seems to come 
from the multigrid algorithm we use to solve for the ini- 
tial data, where linear interpolation is used to prolong 
from the coarse to fine meshes. Linear interpolation in- 
troduces high-frequency components in the fine grid so- 
lution, which is smoothed by relaxation, however relax- 
ation is only applied at interior points. Presumably some 
form of explicit dissipation at parent child boundaries, 
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or higher order interpolation could cure this problem, 
though we find that the dissipation we use during sub- 
sequent evolution is also quite effective at reducing the 
magnitude of this noise. At late times, several grid func- 
tions seem to show anomalously high convergence factors. 
This seems to be due to the fact that the simulations (in 
particular the coarsest resolution one) are not yet that 
close to the convergent regime. To test this would re- 
quire a higher resolution simulation, which would be im- 
practical because of our computer resource limitations 9 . 
However, by looking at an independent residual of the 
Einstein equations, as described next, we can already see 
the trend towards second order convergence using only 
three simulations. 

To check that we are solving the Einstein equations we 
compute an independent residual TZ a p of l[T]l 



= Rap — 47T (2irT a p — g a pT) . 



(62) 



After discretizing the ten residuals TZ a p using the finite 
difference stencils described in the preceding section, we 
compute the residual grid function 1Z at each grid point 
as the infinity norm over the ten residuals. Note that we 
compute (|62f) without reference to the source functions, 
using only the compactified metric elements and scalar 
field. Since we know that 1Z should converge to zero in 
the limit, it is sufficient to compute its convergence factor 
using two resolutions, for example 



2.5 



C? 2 



1.5 



1 

/■ 


II 


ill 


u 

/■' 

7 / 




-■■ 

•»...._._. 


ml 






- /; 

~ a 

II 1 1 III 


ii 


■ Q 

iii 



0.1 0.2 0.3 0.4 0.5 
t 



FIG. 3: Convergence factors l|61|l for representative grid func- 
tions from the simulation described in Sec. IV Al The points 
denote the times when Q was calculated, and correspond to 
times when the entire grid hierarchy was in sync. Note that 
we only show <2§ xx at t = 0, as all the other functions are 
exactly known then, and hence Q is ill-defined. This plot 
shows that the solution is close to second order convergent, 
with some caveats discussed in the text. 



Qn = — {\n\\n 2h \\-\n\\n h \\). 
hi 2 



(63) 



Fig. 0] shows Qk computed using both \R-ih 1 'R-h\ and 
[TZ4h,Ti-2h}- This plot shows that we are tending towards 
second order convergence as the resolution is increased. 



B. Schwarzschild Black Hole Evolution 

In this section we briefly show how well the current 
code can evolve a Schwarzschild black hole in Painleve- 
Gullstrand coordinates. The analytic solution is used 
for initial conditions as described in Sec. IIVFI with 
M = 0.05, and using (|15fl to keep the source functions 
frozen-in during evolution. A ( "lego" ) spherical excision 
region of radius 1.2M was used. We ran three axisym- 
metric simulations, each with identical grid hierarchy, 
though successively higher resolution as described in the 
previous section. The lowest resolution simulation had a 
base grid of 33x17 (spanning — 1..1 in x and 0..1 in y), 
using 6 additional levels of 2:1 refinement, and a Courant 



9 Alternatively, we can choose initial data that is better resolved 
on the coarsest grid; however, the kind of resolution we have here 
is more representative of the resolution we will be able to achieve 
in the near future with the computer power we have access to, 
and so we think this is a fair test of the code. 
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FIG. 4: Convergence factor of the independent residual 
I|62l63fl of the Einstein equations from the simulation de- 
scribed in Sec. IV Al The points denote the times when Q-r. 
was calculated, and corresponds to times when the entire grid 
hierarchy is in sync after an evolution time step (hence there 
are no points at t = 0). This plot shows we are tending toward 
a solution that is second order convergent. 
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factor of 0.125 (the mesh structure is very similar to that 
depicted in Fig[21 however in this example the refine- 
ment is constant in time). To compare, we also ran the 
two lowest resolution simulations of the equivalent prob- 
lem in full 3D; lack of computational resources prevented 
us from running the highest resolution simulation in 3D, 
and, for the same reason, we were not able to run the 
medium resolution simulation as long as the axisymmet- 
ric one 10 . 

As a measure of the accuracy of the simulation, we 
calculate the mass M of the black hole from the area A 
of the apparent horizon: 



M 



A 

16tt' 



(64) 



The mass for the five simulations is shown in Fig. [S] 
Note that we have not calculated any errors associated 
with the numerical integration of apparent horizon area; 
we used the same resolution sphere (33 points in 9, rang- 
ing from to tt) in all cases, hence the error in the area 
calculation will be roughly the same for each run. There 
are a couple of significant things to note from this figure. 
First, even though we were not able to fully compare the 
axisymmetric results with a 3D code, what the partial 
comparison suggests is that explicitly enforcing axisym- 
metry in this case does not have a significant effect on 
the accuracy or runtime of the simulation. Second, even 
though the length of time that we can simulate a black 
hole to within a given accuracy with this code is not too 
long compared to the state of the art these days, the 
trend in increased run-time with resolution is promising. 
In particular, there is not much evidence of exponential 
growth of error early on (though once the error has grown 
to a certain magnitude, the code crashes quickly); rather, 
these plots suggest that the leading order truncation er- 
ror term has somewhere between linear and quadratic 
dependence on time. To see this, let us compare the pu- 
tative runtime of a simulation where the leading order 
error grows exponentially with time, versus polynomial 
growth of the error. For the exponential case, assume 
the norm of the error E(t) as a function of time takes the 
following form 



E(t) = Ch 2 



(65) 



where C is some constant, A is the continuum growth 
factor, and h is the mesh spacing. In other words, this 
situation describes an exponential "constraint violating 



10 Specifically, the medium resolution (2/i) simulation in 3D took 
160 hours of runtime on 128 nodes of the Westgrid Xeon cluster to 
reach t = 55M, using about 120MB of memory on each node. By 
comparison, the highest resolution (h) axisymmetric simulation 
took approximately 240 hours on 24 nodes of UBC's vn4 Xeon 
cluster to reach t = 220Af. In 3D (2D), doubling the resolution 
typically requires 16(8) times the runtime, and 8(4) times the 
memory to evolve to a given physical time in the simulation. 



mode" driven by truncation error terms. Let us solve for 
the evolution time th, to reach a specified error E(t) = E a 
with mesh spacing ft: 



t h 



InEo-lnC -2 In ft 
A 



(66) 



Now consider the following quantity C computed using 
three simulations with differing resolutions: 



th — tlh 
tih — t 4h 



(67) 



Evaluating £ for the case of exponential growth using 
gives 



Ca = 1- 



(68) 



Repeating the calculation for the case of polynomial error 
growth of the form 



E(t) = Ch 2 t p 



gives 







2 2/f> _ 1 

1-2- 2 /p' 



(69) 



(70) 



For linear error growth, £p=i = 4, for quadratic growth 
Cp=2 = 2, and £ p — > 1 in the limit as p — > oo. If we eval- 
uate ( by defining the error to be that in M/Mq from 
FigJSJ using a value of 3% for E we compute £ rss 2.7, 
suggesting polynomial rather than exponential growth. 
However, this number changes as Eq changes (for exam- 
ple setting Eq to 10% suggests faster than exponential 
growth), so we cannot conclusively rule out exponential 
growth. Regardless, from the practical point of view of 
using the current code to investigate black hole physics 
in 3D, we need prohibitively high resolutions to get to 
a useful runtime range of several hundred M, so signifi- 
cantly more work needs to be done to improve the code 
for black hole simulations. 



C. Black Hole Formation 

The final test presented here is gravitational collapse 
of scalar field initial data to a black hole, in axisymmetry. 
To compare with the vacuum black hole simulation of the 
previous section, we used an identical grid structure, and 
chose initial data so that a black hole of roughly the same 
mass (0.05) forms. Specifically, we used a spherically 
symmetric Gaussian pulse (|47|l with 



A 1 = 0.35, A 1 = 0.055, 



(71) 



with the rest of the initial data parameters for $ set to 
zero. We used the same gauge conditions for as de- 
scribed in Sec. IV Al for the convergence test, except here 
the corresponding parameters were kq = 40, £o = 30, 
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FIG. 5: Normalized mass 164H for the evolution of a vacuum 
Mo — 0.05 Schwarzschild black hole in Painleve-Gullstrand 
coordinates. The curve labeled 4/i corresponds to the lowest 
resolution axisymmetrix simulation, while the 2h (h) curves 
are from axisymmetric simulations with twice (four times) the 
resolution. The curves 4/i, 3D and 2h, 3D are from runs with 
identical resolution to the Ah and 2h axisymmetric simulations 
respectively, though the simulations were in full 3D. Note that 
the 2h, 3D simulation curve only extends till roughly t/Mo = 
55 (as we ran out of computer time then), and is effectively 
hidden behind the other curves as M/Mo ~ 1 up till then. 

n = 5, ti = 1/80 and oq = 1. Note that the results 
are not very sensitive to this particular choice of gauge 
parameters; the rule of thumb is that kq and £o of order 
1/A , t\ of order A 1 and n of order unity works reason- 
ably well. Fig. shows the corresponding plot of appar- 
ent horizon mass versus time. The black hole forms af- 
ter about 2M of evolution, after which some accretion of 
scalar field occurs, causing the mass to grow by a bit early 
on. Note also that once we detect an apparent horizon, 
we excise a spherical region 60% the size of the horizon, 
so approximately at 1.2M , again for comparison with the 
previous evolution. At a given resolution this simulation 
(as judged by the mass estimate) has less accuracy com- 
pared to the corresponding vacuum simulation, however 
the trend of increased accuracy with increased resolution 
is roughly the same. 



VI. CONCLUSION 

We have described a new computational scheme for 
numerically solving the Einstein equations based upon 
generalized harmonic coordinates. This extends earlier 
work of Garfinkle |20| . and in some respects is similar 
to the direction pursued by Szilagyi and Winicour [2l| . 



FIG. 6: Normalized mass I64H from the axisymmetric evolu- 
tion of a black hole formed via the gravitational collapse of a 
scalar field. The value of Mo used was the largest, convergent 
value of M from the three simulations, which is a reasonable 
estimate of the final mass of the black hole. For comparison, 
the grid structure used for the simulations was identical to 
the corresponding axisymmetric simulations shown in Fig|^] 



Some of the topics covered included suggestions for im- 
posing dynamical gauge conditions, a new technique of 
implementing the Cartoon method 32] for simulating ax- 
isymmetric spacetimes with a Cartesian code, a direct 
discretization scheme for second-order-in-space-and-timc 
partial differential equations, and the use of a spatially 
compactified coordinate system. One attractive feature 
of harmonic evolution is that the principal part of the 
Einstein equations reduce to wave equations for each met- 
ric element. This, together with the use of a second or- 
der discretization scheme, keeps the number of variables 
and constraint equations to a minimum, and the hope is 
that this will make it easier to achieve stable evolution. 
The use of a spatially compactified domain allows one to 
impose correct asymptotic boundary conditions for the 
simulation, and thus we automatically have constraint 
preserving boundary conditions. The advantage of our 
Cartoon method over the original is that no interpola- 
tion is needed, and the simulation is performed on a 2D 
slice of the spacetime, thus simplifying the process of in- 
corporating the code into an adaptive mesh refinement 
framework. Furthermore, the technique is not particular 
to finite difference codes, and can be used with spectral 
methods, for instance. 

Preliminary test simulations of black hole spacetimes 
suggest that this scheme holds promise for being appli- 
cable to many problems of interest, including the binary 
black hole problem, black hole-matter interactions, and 
critical gravitational collapse. However, a lot of research 



15 



still needs to be done, at both the analytical and numer- 
ical levels, before this scheme may produce new physical 
results. In particular, it would be useful to analyze the 
mathematical well-posedness of the fully discrete system, 
including a variety of possible gauge evolution equations. 
The majority of techniques for analyzing hyperbolic sys- 
tems require reduction to first order form (recently sim- 
ilar techniques have been develo ped for second order in 
space, first order in time systems |58l \E$L l60| : also, in [6l| 
the BSSN system is analyzed by converting to first order 
form, however the constraints introduced by this reduc- 
tion are shown to obey a closed evolution system that is 
independent of the other constraints, implying that the 
original second order system is well-posed). At the nu- 
merical level, a broader class of initial conditions needs 
to be explored, such as black hole/matter interactions 
and black hole collisions. This is of course one of the 
primary long term goals of the code, however early tests 
indicate that a significant number of adjustments and im- 
provements (to dissipation and extrapolation operators, 
for example) may be needed, in addition to more sophis- 
ticated gauge conditions than discussed here, before such 
scenarios could be simulated with sufficient accuracy and 
length of time for new results to be obtained. 
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APPENDIX A: EVOLUTION OF 
AXISYMMETRIC SPACETIMES 

As described in Sec. II V CI we can efficiently simu- 
late axisymmetric spacetimes along a single z — slice 
of the spacetime by replacing all z gradients in the field 
and matter evolution equations with appropriate x and y 
gradients, as dictated by 153|) . Here we list the equations 
for gradients of the regular components of the metric 
and scalar field $ with respect to the compactified coor- 
dinates (Sec. 



IIV Bll . and give the corresponding on-axis 
regularity conditions. 

The first z derivatives are: 



d z gxy\z=o = -g~xz(n/2y) 

d z 9xz\z=o = g~xy(ir/2y) 

dzg yy \z=o = -3yz(n/y) 

d z gyz\z=o = (g yy - gzz)(Tr/2y) 

d z g~zz\z=o = 9yz(^/y) 



(Al) 



Mixed z — t, z — x and z — y second derivatives are calcu- 
lated by taking the appropriate derivative of (jAljl . Sec- 
ond derivatives with respect to z are computed as follows: 



d z d z g. 



'af3\z=0 



( d y 9afi I^Caf. 



2\y(l + f) 2f 



d z d z <P\ z =o 



2y(l + y 2 ) 
where the coefficients C a p are 

C tt = 

Ct S = o 

C ty = —gt y 

Ctz = —gtz 

Cxx = 
Gxy — Qxy 
Cxz 9xz 

Cyy = 2(§z Z — (jyy) 

C y z = —^g y z 

Cgg = 2((jyy — (] ZZ ) 

The on axis regularity conditions are 

d y g t t\ y =o = 
d y g~tx I y—o = 

9ty\y=0 = 



(A2) 



(A3) 



9ts\ y =o 







d y 9xx\y=t) — 

9x y \y=0 = 

gx Z \ y =o = 

dy9yy\v=o = 

9y9yz\y=0 = 

d y g zz \ y =o = 

fy*lv=o = 

9yy\y=o = 9zz\ y =o 



(A4) 



d z 9tt\z=o 
d z 9tx\z=o 
d z gt y \z=o 
d z gtz\z=o 
d z gxs\z=o 






-g tz {-K/2y) 

g ty {-n/2y) 





To compute the z gradients and regularity conditions 
for the source functions in the code, we simply substitute 
the results from the calculation for the metric into the 
definition of the source functions 17I37|) . 
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APPENDIX B: STABILITY ANALYSIS OF A 
SECOND ORDER IN SPACE AND TIME 
EVOLUTION SCHEME 



where 



£ = sin(fcAx/2), 



(B9) 



Here we give a von Neumann-like stability analysis of 
the one dimensional flat space wave equation using the 
discretization scheme described in Sec. IIVI Second or- 
der in time schemes for the wave equation are not very 
common in the literature, so the example given here is 
to demonstrate that the method is inherently stable, ig- 
noring the complications of boundaries, excision, non- 
constant coefficients and non-linear lower order terms of 
the full problem (the analysis of which is beyond the 
scope of this paper). Even though dissipation is not 
needed in this example, we add it as applied in the code 
to demonstrate how it works. 

The model wave equation for <&(x, t) is 

*,« - = 0. (Bl) 

Discretization of this equation using the stencils in Tab. 
ID gives 

$«+i_2$«+$™- 1 -A 2 - 2$J + = 0, (B2) 

where $" = $(x = jAx,t — nAt) and A = At/ Ax is 
the Courant factor. This immediately gives an explicit 
update scheme for the unknown <E> ra+1 given information 
at two past time levels, 



$ n+1 = 2$" - $ 

3 3 3 



a 2 - + $?_ x ) (B3) 



(note that the iterative relaxation method described in 
Sec. IIVI gives exactly the same update scheme in this 
case). It is mathematically simpler to analyze this equa- 
tion using an equivalent two time level scheme by intro- 
ducing the variable 



viz™ = A"" 1 

3 3 ' 



after which (IB3I) becomes 



$«+! = 2$" — + A 2 ($"+1 — 2$™ + 



3 



<J>" 
3 ' 



(B4) 



(B5) 



As (lB5jl is linear with constant coefficients, we can com- 
pletely characterize its stability properties by analyzing 
the evolution of individual Fourier modes of the form 
c{t)e lkx . To this end, let 



$(x, t) ee a{t)e lkx 
*(x,t) ee b(t)e tkx . 



Substituting this into (|B5(I gives 



I = Za — b — AA t, a 

un+l _ n n 



(B6) 
(B7) 



(B8) 



and we have used the identity —4 sin 2 (fc Ax/2) = 
e -ikAx _ 2 + e lkAx . Note that the smallest wavelength 
that can be represented on a numerical grid is 2Ax (the 
Nyquist limit), which corresponds to a largest possible 
wave number k — tt/Ax, hence £ ranges from to 1. 

As described in Sec. IIV A II we apply numerical dissi- 
pation to all past time level variables, prior to the update 
step, by first calculating the high-frequency component 
of the function using (|30|l . and then subtracting it from 
the function via (|31|l . For a grid function /" 



the high-frequency component rj™ is defined as 



1 

16 

m 4 , 



(/jL 2 -4/jL 1 +6/;-4/JV 1 + /JV 2 ) 



(BIO) 



and filtering amounts to modifying /" as follows 



f n 

J 3 



= /;(i~^ 4 ) 



(Bll) 



where e ee 1 - e^ 4 . As C 6 [0..1] and e e [0..1], e G [0..1]. 
With this form of dissipation (which is linear, and hence 
fits into the Fourier analysis of the evolution scheme) 
applied to both <I>™ and 1B8II becomes 



,n+i _ 



a-' = e [2a™ (l - 2A 2 £ 2 ) - b n ] 
b n+1 = ea n . 

In matrix form, the update step can be written as 



(B12) 





n+l 




n 


a 


= A 


a 




b 




b 





where 



A = 



2(1-2A 2 C 2 ) -1 
1 



(B13) 



(B14) 



The numerical evolution will be stable if the eigenvalues 
A± of A all lie on or within the unit circle in the complex 
plain. A straight-forward calculation gives 



1 - 2£ 2 A 2 ± ^Av 7 ! - £ 2 A 2 



(B15) 



The expression within the square root of (|B15|) is strictly 
non-negative if we require that A € [0..1]. The magnitude 
of the eigenvalues are 

l|A±||=6. (B16) 
Hence, for e < 1 the numerical scheme is stable; in fact, 
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without dissipation (e = 1) the scheme is inherently sta- ble and non-dissipative. 
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